A Gentle introduction to Markov Process

A Hands on Python Codes in the case of Inventory Managements

Intro

If you are an analytics professional (whether a Data Scientist, Analyst, etc) you will be dealing often with Processes where the process is time indexed and follows an uncertain path. Imagine, for example, that you are a data scientist working at an Energy company. Your task involves keeping track of uncertain path of commodity price.

The commodity price (Oil for instance) will follow a path at time steps $t=0, 1, 2, \cdots$. This is a uncertain Process where it was indexed by time. However, if we want do analysis and make decision you need to internally represent the process.

Think of State as tool to internally represent the uncertain process. Lets back to example of Oil Price. Let's say today's oil price is $100. I can represent that information with saying $S_0= 100\$$, then tomorrow the price will be different, can be represented as $S_1$ and this sequence continues $S0, S_1, S_2,...$.

Here in this short log, we want to talk anbout Markov Process. As you can from the name, Markowv Process deal with a uncertain path. It mean Markov Process is sequence of random states $S_t \in \mathcal{S}$, as at time $t=0,1,2,3,\cdots$, we can reprsent the process with:

$$S_0, S_1,S_2,\cdots,S_{t-1},S_{t}$$

Markov Property

We call the process has Markov Property where the transition of states has following feature:

$$\mathbb{P}[S_{t+1}|S_t,S_{t-1},\cdots, S_0] = \mathbb{P}[S_{t+1}|S_t]$$

Let me explain the above equation in simple term. Lets say you are playing game that you move between three spaces, A, B, C. And Let's say you are in space A, there is 50% chance you will stay at A in next stop, 20% you will move to state B and 30% you move to state C. This means that , it does not matter what was history of your past moves, as long as you are at space A ($S_t= A$), the probality of next spaces $P(S_{t+1}=A) = 0.5$, $P(S_{t+1}=B) = 0.2$ and $P(S_{t+1}=C) = 0.3%$.

Note that in above example, we only care about the transition probability of process, represented with State. Here, we do not care about the decision or reward (As in this blog I am only talking about Markov Process). I see this mistake a lot, Markov Process is confused with Markov Reward Process or Markov Decision Process.

Now, let's do some coding to understand what a Markov Process is. Here, I'll give an example of managing the inventory of a bike shop, and we'll see how changes in the store's inventory can be internally represented as "State".

A Basic Example of Bike Shop Inventory Management

display(Image(filename="img/bikeshop.jpg"))

#reference : 
#https://unsplash.com/photos/zbUH21c9ARk?utm_source=unsplash&utm_medium=referral&utm_content=creditShareLink

Imagine you own a bike shop with an inventory that can only hold a certain number of bicycles (limited capacity). Every day, some customers come to buy bikes from your shop. For example, Let's assume your shop can hold a maximum of 5 bikes. Assume that today (Wednesday) and you count and find that you have three bikes in your shop. However, you also know that tomorrow (Thursday) there will be some demand, meaning some of your bikes will be sold. The exact demand for bikes tomorrow is uncertain, and you're not sure how many costumers you will have.The challenge is you need to figure out how many bikes to order each day. You shouldn't order too many because keeping a lot of bikes in your inventory costs money. On the other hand, you must order enough to meet the potential demand of your bike shop, (if not you have missed the demand, bad for your business).

Markov Process

Example of Bike Shop Inventory Management

In this blog, our main focus is on the Markov Process. It means in this blog, wil have (assume) fixed policy(which, in this case, policy is how much bike to order every day). To model the Markov Process of this example, first we need to frame this problem (what is state?) , secondly we need to build the model of state transition.

In this problem, the state (the internal representation of process) can be described by two variable:

• $\alpha$ : is number of bike you already have in store (On-Hand Inventory at 6pm)
• $\beta$ : is the inventory on a truck from the supplier (that you had ordered teh previous day) and will arrive next day.

The sequence of events in 24-hr cycle for this bike shop is as follow:

• $\alpha$ : is the inventory in the store (On-Hand Inventory at 6pm)
• $\beta$ : is the inventory on a truck from the supplier (that you had ordered teh previous day)
• Observe State $S_t: (\alpha, \beta)$ at 6pm store closing
• Order Quantity: = $max(C-(\alpha + \beta), 0)$ immdeiately after 6pm
• Receive inventory at 6am if you had ordered 36 hrs ago
• Open the store at 8pm
• Experience random demand $i$ with poisson probabilities: $$f(i) = \frac{e^{-\lambda}\lambda^i}{i!}$$
• Close the store at 6pm and observe new state $S_{t+1}: (\alpha, \beta)$

Similar to the explanation above, the $S_t = (\alpha_t + \beta_t)$ is internal representation of this stochastic process, and in Markov Process, we want to model how this stochastic process evolves. Let me give example. Let's Say David is owner of Bike shop. On a Wednesday (6PM), he has 2 bikes in his store, and 1 bike he ordered on Monday 6PM. His state would be:

$$S_{t= \text{Wednesday}} : (\text{2 bike in shop at Wednesday 6PM} + 1 \text{bike that on truck will arrive tomorrow (Thursday 6AM}) = (\alpha = 2, \beta = 1) $$

So the state at the time Wednesday 8 PM is $S=(\alpha=2, \beta = 1)$. Then, each day there is a random (non-negative integer) demand for bikes, with demand modeled following Poisson distributions (with Poisson parameter $\lambda \in \mathbb{R}_{>0}$). A demand of $i$ bikes for each $i=0,1,2\cdots$ occurs with probability:

$$\mathbb{P} (\text{experiencing a demand of i on a given day})=f(i) = \frac{e^{-\lambda}\lambda^i}{i!}$$

We can visualize this distribution to see probability of experiencing different demands, given we pick the $\lambda =1$

#from rich import print, pretty
#pretty.install()

from rich import pretty
pretty.install()
from IPython.display import Image, display


# import matplotlib and some desired styling
import matplotlib.pyplot as plt
%matplotlib inline

plt.style.use("ggplot")
plt.rcParams["figure.figsize"] = [10, 6]

# need numpy to do some numeric calculation
import numpy as np

# poisson is used to find pdf of Poisson distribution 
from scipy.stats import poisson

x_values = np.arange(0,10)

# pdf Poisson distri with lambda = 1
pdf_x_values = poisson.pmf(k=x_values,mu=1)

plt.bar(x_values, pdf_x_values, edgecolor='black')
plt.xticks(np.arange(0, 10, step=1))
plt.xlabel("Costumer Demand")
plt.ylabel("Probability of Costumer Demand")
plt.title("Probability Distribution of Costumer Demand with $\lambda = 1$ ")
#plt.savefig("fig/poisson_lambda_1.png", dpi=300)
plt.show()

Constructing The Probability Transition of the Markov Process:

Now that we have a initial understanding of what a state is, $S_t = (\alpha, \beta) $, the uncertainty that introduces the concept of a stochastic process is the customer demand ($i$). We can build upon on how the probability transition evolves within this problem. I will do coding for that one, but first let's explain it in simpler term.

The Probability transition of this problem has two cases, case 1 and 2.

• Case 1)

If the demand $i$ is smaller than the total inventory available on that day, $\text{initial inventory}= \alpha + \beta$

$\mathbb{P}[\text{State} \; (\alpha, \beta) \rightarrow \; \text{State}(\alpha + \beta -i, C-(\alpha, \beta))]$ = $f(i) \quad \quad \quad 0<i<\alpha+\beta -1 $

• Case 2)

If the demand $i$ is more than total inventory availabe at that day, $\text{initial inventory}= \alpha + \beta$

$\mathbb{P}[\text{State} \; (\alpha, \beta) \rightarrow \; \text{State}(0, C-(\alpha, \beta))]$ = $\sum_{j=\alpha + \beta}^{\infin}f(j)=1-F(\alpha+\beta -1) \quad \quad \quad i>\alpha+\beta $

Where $F(\alpha + \beta -1)$ is the CDF of Poisson distribution. Having understanding the background of this problem, now we can write some python code to understand better. The crucial aspect of this python code is designing the data structure to store the transition of Markov Process. To do that I am designing the data structure as a Dictionary I call, "MarkovProcessDict". The keys of this dictionary correspond to the current state, and the values (represented in dictionary) are the next states along with the probabilities associated with each next state. This is sample example how the data structure MarkovProcessDict look like:

from typing import Dict

MarkovProcessDict = {"Current State A":{"Next State 1, from A": "Probability of Next State 1, from A",
                                        "Next State 2, from A": "Probability of Next State 2, from A"},
                    
                     "Current State B":{"Next State 1, from B": "Probability of Next State 1, from B",
                                        "Next State 2, from B": "Probability of Next State 2, from B" }}

MarkovProcessDict

{
    'Current State A': {
        'Next State 1, from A': 'Probability of Next State 1, from A',
        'Next State 2, from A': 'Probability of Next State 2, from A'
    },
    'Current State B': {
        'Next State 1, from B': 'Probability of Next State 1, from B',
        'Next State 2, from B': 'Probability of Next State 2, from B'
    }
}

In the above code 'Current State A' is the initial state, that we can move to two possible states 'Next State 1, from A' and 'Next State 2, from A' , each with probability of Probability of Next State 1, from A and Probability of Next State 2 from A.

MarkovProcessDict: Dict[tuple, Dict[tuple, float]] = {}

user_capacity = 2
user_poisson_lambda = 1.0

# We are condiering all possible states
# That we can face in running this bike shop
for alpha in range(user_capacity+1):                            
                                                               
    for beta in range(user_capacity + 1 - alpha):
        
        # This is St, the current state
        state = (alpha, beta)                                   

        # This is initial inventory, total bike you have at 8AM 
        initial_inventory = alpha + beta                         
        
        # The beta1 is the beta in next state, irrespctive of current state (as the decsion policy is constant)
        beta1 = user_capacity - initial_inventory
        
        # List of all possible demand you can get
        for i in range(initial_inventory +1):

            # if initial demand can meet the demand
            if i <= (initial_inventory-1):
                
                # probality of specifc demand can happen
                transition_prob = poisson.pmf(i,user_poisson_lambda)
                
                # If we already defined teh state in our data (MarkovProcessDict)
                if state in MarkovProcessDict:
                    
                    MarkovProcessDict[state][(initial_inventory - i, beta1)]= transition_prob
                
                else:
                    
                    MarkovProcessDict[state] = {(initial_inventory - i, beta1):transition_prob }
                         
            # if initial demand can not meet the demand
            else:
                
                # probability of not meeting the demands
                transition_prob = 1- poisson.cdf(initial_inventory -1, user_poisson_lambda)

                if state in MarkovProcessDict:
                    
                    MarkovProcessDict[state][(0, beta1)]= transition_prob
                    
                else:

                    MarkovProcessDict[state] = {(0, beta1 ):transition_prob }

MarkovProcessDict

{
    (0, 0): {(0, 2): 1.0},
    (0, 1): {(1, 1): 0.3678794411714424, (0, 1): 0.6321205588285576},
    (0, 2): {
        (2, 0): 0.3678794411714424,
        (1, 0): 0.3678794411714424,
        (0, 0): 0.26424111765711533
    },
    (1, 0): {(1, 1): 0.3678794411714424, (0, 1): 0.6321205588285576},
    (1, 1): {
        (2, 0): 0.3678794411714424,
        (1, 0): 0.3678794411714424,
        (0, 0): 0.26424111765711533
    },
    (2, 0): {
        (2, 0): 0.3678794411714424,
        (1, 0): 0.3678794411714424,
        (0, 0): 0.26424111765711533
    }
}

for (state, value) in MarkovProcessDict.items():
    print("The Current state is: {}".format(state))
    for (next_state, trans_prob) in value.items():
        print("The Next State is  {} with Probability of {:.2f}".format(next_state, trans_prob))

    #print(state)
    #print(value)

The Current state is: (0, 0)
The Next State is  (0, 2) with Probability of 1.00
The Current state is: (0, 1)
The Next State is  (1, 1) with Probability of 0.37
The Next State is  (0, 1) with Probability of 0.63
The Current state is: (0, 2)
The Next State is  (2, 0) with Probability of 0.37
The Next State is  (1, 0) with Probability of 0.37
The Next State is  (0, 0) with Probability of 0.26
The Current state is: (1, 0)
The Next State is  (1, 1) with Probability of 0.37
The Next State is  (0, 1) with Probability of 0.63
The Current state is: (1, 1)
The Next State is  (2, 0) with Probability of 0.37
The Next State is  (1, 0) with Probability of 0.37
The Next State is  (0, 0) with Probability of 0.26
The Current state is: (2, 0)
The Next State is  (2, 0) with Probability of 0.37
The Next State is  (1, 0) with Probability of 0.37
The Next State is  (0, 0) with Probability of 0.26

import graphviz  

d = graphviz.Digraph(node_attr={'color': 'lightblue2', 'style': 'filled'}, 
                    )
d.attr(size='10,10')

d.attr(layout = "circo")

for s, v in MarkovProcessDict.items():

        for s1, p in v.items():
            
            d.edge("(\u03B1={}, \u03B2={})".format(s[0],s[1]), 
                   "(\u03B1={}, \u03B2={})".format(s1[0],s1[1]), label=str(round(p,2)), color="red")
            
print(d)

digraph {
	node [color=lightblue2 style=filled]
	size="10,10"
	layout=circo
	"(α=0, β=0)" -> "(α=0, β=2)" [label=1.0 color=red]
	"(α=0, β=1)" -> "(α=1, β=1)" [label=0.37 color=red]
	"(α=0, β=1)" -> "(α=0, β=1)" [label=0.63 color=red]
	"(α=0, β=2)" -> "(α=2, β=0)" [label=0.37 color=red]
	"(α=0, β=2)" -> "(α=1, β=0)" [label=0.37 color=red]
	"(α=0, β=2)" -> "(α=0, β=0)" [label=0.26 color=red]
	"(α=1, β=0)" -> "(α=1, β=1)" [label=0.37 color=red]
	"(α=1, β=0)" -> "(α=0, β=1)" [label=0.63 color=red]
	"(α=1, β=1)" -> "(α=2, β=0)" [label=0.37 color=red]
	"(α=1, β=1)" -> "(α=1, β=0)" [label=0.37 color=red]
	"(α=1, β=1)" -> "(α=0, β=0)" [label=0.26 color=red]
	"(α=2, β=0)" -> "(α=2, β=0)" [label=0.37 color=red]
	"(α=2, β=0)" -> "(α=1, β=0)" [label=0.37 color=red]
	"(α=2, β=0)" -> "(α=0, β=0)" [label=0.26 color=red]
}

d

d.render('img/MP', format="pdf")

'img/MP.pdf'